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Spatial correlation functions provide a glimpse into the quantum correlations 
within a quantum system. Ions in a linear trap collectively form a nonuniform, 
QQ discretized background on which a scalar field of phonons propagates. Trapped ions 

^^ have the experimental advantage of each having their own "built-in" motional de- 

^ tector: electronic states that can be coupled, via an external laser, to the ion's 

,__i vibrational motion. The post-interaction electronic state can be read out with high 

1^ efficiency, giving a stochastic measurement record whose classical correlations reflect 

pg the quantum correlations of the ions' collective vibrational state. Here we calcu- 

CN late this general result, then we discuss the long detection-time limit and specialize 

,— , to Gaussian states, and finally we compare the results for thermal versus squeezed 

i-^ states. 

^ I. INTRODUCTION 

I In contexts as diverse as cosmology, quantum optics and condensed matter physics, spa- 

>■ tial correlation functions provide experimental access to important features of spatially dis- 

^ tributed quantum systems. In a cosmological setting, spatial correlations in the temperature 

ly-^ distribution of the cosmic microwave background provide direct access to the fluctuations 

cn of a primordial quantum field during infiation. Details of these spatial correlations provide 

t^^ direct and detailed tests of infiationary cosmological models [T]. For example, Grishchuk [2] 
has suggested a model in which anisotropies refiect the underlying statistics of squeezed 
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(^ vacuum states. In condensed matter physics, the change in the length scale of correlations 



> 
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of the XY model in two dimensions can reveal a Kosterlitz-Thouless transition and vortex 
phases [Hll]. The spatial correlation functions of a Bose-Einstein condensate, as reveled by 
light scattering from a freely expanding condensate reveal details of the quantum state of 
^ the condensate before expansion [H]. As one example of this, mean field energy of a con- 

densate is a direct measure of the second order correlation function [Hj. The atom loss-rate 
due to three-body recombination in a BEG is directly related to the probability of finding 
three atoms close to each other and can therefore act as a probe of a third-order correlation 
function [7]. In this paper we show that similar two-point spatial correlation functions can 
be used to probe the collective vibrational state for a string of trapped ions [8] . Our work 
is related to that of Franke- Arnold \9\, which considers the spatial coherence properties of 
just two harmonically trapped particles. 

Using external lasers, it is possible to couple an internal electronic transition to the 
vibrational motion of the ion [TU]. Indeed this is how laser cooling of multiple ions to the 
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vibrational ground state of one or more of their collective normal modes of vibration |TT] . 
Simple quantum information processing tasks can be realised by coupling internal states 
through collective motional degrees of freedom fTIi IT3] . In this paper, we show how the 
ability to couple internal and vibrational degrees of freedom, plus the ability to efficiently 
read out the internal electronic state using a fluorescence shelving scheme [S], provides a 
path to experimentally measuring spatial correlation functions for the collective vibrational 
degrees of freedom. This provides a discrete analogue of spatial correlation functions for a 
scalar quantum field. 

The scheme has two components: (1) external lasers weakly couple the motion of two 
distinct ions to their local displacement from equilibrium; subsequently, (2) external readout 
lasers are used to probe the electronic state of all the ions. In the second stage, some ions 
will "light up" (made more precise below), while others remain dark. Thus, for N ions, 
we get an iV-element stochastic binary string (ei,e2, . . . ,6^), with e^ = 1 if an ion lights 
up and Cm = otherwise. This is the measurement record, with the ordering of the string 
corresponding exactly with the position ordering of the ions in the trap. We can then define 
the empirical two-point correlation function Pmn as the classical average Pmn = £[^m^n\- 
Our objective in this paper is to relate this classical correlation function to the quantum 
mechanical correlations of the ion displacement amplitudes (e.g., terms such as {qm<in)) 
and, furthermore, to determine characteristic experimental signatures of different collective 
quantum states of motion. 

II. SPATIAL CORRELATION FUNCTIONS 

We summarise the standard results for the spatial correlation functions of a scalar field 
(see, for example, Ref. [H]). A scalar field 0(x, t) may be expanded as 

0(x, t) = ^ aktik(x, t) + aj^Mk(x, t) , (2.1) 

k 

where Mk(x, t) is a positive-frequency mode function, and mJ^(x, t) is its negative-frequency 
counterpart. In the commonly used case of plane waves (with box normalization), we can 
define the positive- and negative-frequency components of 0(x, t) as 

V;(x,t):=-^5^ake+^('^— '=*) (2.2a) 

and 7/^t(x,t) := ^^4e-^('^-^-^'=*) , (2.2b) 

k 

respectively. In the continuum limit, these operators satisfy the standard boson equal-time 
commutation relations: 

[V^(x,t),V^t(x',t)] = 5(x-x'). (2.3) 

With all of the time-dependence in the exponential e^*'^'=*, we can change to the Schrodinger 
picture simply by simply removing this piece (which is equivalent to evaluation at t = 0). 

Now in the Schrodinger picture, we define the first- and second-order normally ordered 
spatial correlation functions, 

G(i)(x,x') :=(^^(x)^(x')> (2.4a) 

and G(2)(x,x') := (^Hx)^t(x')^(xXx)) . (2.4b) 



We also define the normalized correlation functions, 
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If each mode is independently excited to a thermal state, with no mode-mode correlations, 
we can show that [H] 

^(2)(x,x') = l + |^«(x,x')f , (2.6) 

where the second term represents bosonic bunching. In the continuum limit. 



^(i)(x,x') = ^ f d^k n(k)e''^-("-"') 



(2.7) 



where ?T,(k) represents the thermal occupation number of mode k and N is the total occu- 
pation number, given by 

N:= f d^knik) . (2.8) 

It is apparent that (^'■^•'(x, x'), as a function of (x — x'), starts at a value of 2 and decreases 
to 1. How fast it decreases is a measure of the range of spatial correlation functions and 
depends on the k-dependance of n(k). 

The example of the thermal states is a special case of a general result for all classical 
states of the field for which 

(?(2)(x,x)>(7(2)(x,x + y) (2.9) 

Nonclassical states violate this inequality. The particular example of squeezed light has 
been studied in some detail [15] . The effect of the squeezing is to induce spatial correlations 
that modulate an effective thermal background density. In some cases this reduces the 
correlation function G^'^\'x, x') below the thermal value of 2, a phenomenon related to photon 
antibunching in the optical case ^2]- It can also increase it above 2. The ability to change 
the first and second order correlation functions in this way is what lies behind the new field 



of quantum imaging [T7|. In Section VI, we will contrast thermal and squeezed states in the 
case of vibrations in a linear ion trap and show that the latter have a similar nonclassical 
signature. 

III. MEASUREMENT OF ION TRAP SPATIAL CORRELATIONS 

A. Normal modes of vibration 

Linearizing the potential created by the overall harmonic potential due to the trap elec- 
trodes and the mutual Coulomb repulsion between the ions, we can model the collective 
motion of the ions in a linear trap as a collection of coupled harmonic oscillators. With a 
coupling matrix A obtained from the linearized combination of these potentials, the Hamil- 
tonian is given by [TT] 

1 ^ Mi/2 



where q = (gi, . . . , QnY' is a column vector of operators corresponding to the displacement of 
each ion from its equilibrium position, p = (pi, . . . ,Pn)'^ is a column vector of corresponding 
momenta, M is the mass of each ion, and u is the effective harmonic trap frequency provided 
by the trap electrodes (typically, i/ ~ a few MHz [5]). The coupling matrix A is symmetric 
and positive-definite. Thus, we can diagonalize it: 

A = B^AB , (3.2) 

where A = diag(/ii, . . . , ^in) is a diagonal matrix of positive eigenvalues, arranged in as- 
cending order. The orthogonal matrix B defines the transformation to normal-mode coordi- 



nates Q = Bq and momenta P = Bp. Equation (3.1) can be rewritten in these coordinates 
as 

1 rr Mz/2 ^ 

^2M^ 2 ^P 
p=i 

= 22^'^p^l^p ' (3-3) 

p 

where 

Up := u^/JI^ (3.4) 



is the oscillation frequency of normal mode p, and (since the normal modes oscillate in- 
dependently) we have diagonalized this Hamiltonian using the standard prescription for a 
set of independent harmonic oscillators, ignoring the zero-point energy. The normal-mode 
raising and lowering operators satisfy 

[ttp, a],,] = Spp> , (3.5) 

as is appropriate for independent bosonic modes. 

The local oscillations qm about the ions' equilibrium positions are given, in the interaction 
picture, by pjj 




(3.6) 

where m E {1, . . . , N} labels the ion, and bm is an entry of B, which defines the spatial mode 
functions of the normal mode.^ We can define a unitless version of this local displacement 
oeprator, as well: 

^ l(p) 

<l^mit) := E -^.(^'"""'^P + ^'"^-'^D ■ (3-7) 

p=i ^P 



^ Prefactor and phase conventions for qm vary in the hterature. 



The connection between the two is given by 



^m{t) = ^ ^<Pm{t) ■ (3.8) 

The positive- and negative-frequency components of (f)m{t) niay be defined analogously to 



those of Eqs. (2.2): 



p=i ^^p 



'"-'ar, 



^ iXp) 
and ^Ut) := ^ ^e+^^^aj . (3.9) 

p=i ^P 

Because the coupling matrix A only approximates a simple "balls on a string" coupling 
(characteristic of a bosonic field), the equal-time commutation relations for these operators 



only approximate those of a boson field (compare with Eq. (2.3)): 



^ hip) hip) 



l^UtUUt)] = E -^ = {^-"'")n.n ■ (3-10) 



p=i \/^^p 



Still, the fact that the ions' normal-mode operators commute properly, as in Eq. (3.5), means 
that we don't need to worry about this, as long as we do our calculations using the normal 
modes explicitly. 



B. Laser-induced coupling of vibrational and electronic states 

The abstract relations between the field spatial correlation functions described in Sec- 
tion [IT] ultimately are manifest in the observed correlations in spatially distributed detectors 
of some kind. In the case of quantum optics, for example, one imagines that the field falls on 
a photodetector array. Each element of the array produces a photocurrent I{x,t) indexed 
by the position of that photodetector in the array. We can then look at cross-correlations 
between photo currents from different detectors in the array. It is the the objective of 
measurement theory to determine how those spatially dependent current-current correlation 
functions are determined at a fundamental level by the field spatial correlation functions 
themselves. Explicit formulae are given in Kolobov [T3] using the standard quantum optics 
theory of photodetection. 

Our measurement model is different from that considered in quantum optics and so 
we now explicitly make the connection between the observations and the underlying field 
correlation functions for the motion of the ions in the trap. The feature of our model is 
that the internal electronic state of each ion can be turned into a local detector for the 
displacement of that ion. To achieve this for a given ion, its internal state must become 
correlated with its linear displacement from equilibrium. 

This correlation is provided by an external laser, which is used to drive an elec- 
tronic transition between two meta-stable electronic levels \g) and |e), separated in energy 
by hujji^ fHl [m HB] . The interaction between an external classical laser field and the mth ion is 



described, in the dipole and rotating-wave approximation, by the interaction-picture Hamil- 
tonian [H [11] 



r(m) 



Am) (,\ik COS eqm{t) _ Jm)^,-. ik cos eqm{t) 



Hp = -thno a'^'{t)e'''°'"^^^'^ - ar'{t)e-""'°'"'"^^'^ , (3.11) 

where Qq is the Rabi frequency for the laser-atom interaction (typically, Qo ~ 100 kHz [8]), 
k is the magnitude of the laser's wave vector k, which makes an angle 6 with the trap 
axis, qm{t) is the interaction-picture position operator for the mth ion, and cr^ (t) are its 
interaction-picture electronic raising and lowering operators. Explicitly, 

af ^t) = e±^^Vf\ (3.12) 

where a^"^ = \e)^Jg\ and a'l'^ = |^)„Je|, and 

A = ti^A — cul (3.13) 

is the detuning of the laser below the atomic transition. The size of the rms fluctuation 
in Qm as compared to the wavelength of the laser is measured by the Lamb-Dicke parameter 



where ^Jh|2Mv is the rms fluctuation of the center-of-mass mode of the ions in the ground 
state. This quantity is representative of the overall rms fluctuations of the mth ion, de- 
noted Axrms, since the frequencies v^ for all higher normal modes of the ions remain within 
an order of magnitude of v for small numbers of ions (up to A^ ~ 10), realistic for cur- 
rent experiments [TT]. Typical values of the Lamb-Dicke parameter are r/ ~ 0.01 to 0.1 [19] . 
When 7] -C 1, the so-called "Lamb-Dicke limit," the ion is well localized with respect to the 



wavelength of the laser, and we can expand the exponentials in Eq. (3.11 ) to flrst order in 77, 
which is equivalent to flrst order in k cos 9 qm(t), giving 

H\"'\t) ~ hnoa^'^\t) + mokcos9qmit)a^r\il , (3-15) 

carrier sideband 

where a^^(t) = e+^^Vf )+e-*^Vi"\ and4'"^(t) = -te+^^'a^^^ +te-'^'a^^\ The flrst term 
corresponds to excitation of the transition directly by the laser, while the second couples 
the atomic transition to vibrational motion. 

The "carrier" term is resonant when A = and corresponds to direct excitation of the 
atomic transition, which does not couple to the vibrational motion at all. For sufficiently 
long-time detection, another rotating-wave approximation may be made in which only res- 
onant (nonoscillatory) terms are kept. In this case, if the carrier transition is sufficiently 
off- resonant (A 7^ 0), then it can be neglected, leaving 

H^r\t) ^^^^^ ^^or] y ^Ae-'^^'ar + e^^^*at)(e^^Vi") + e'^^V^'")) . (3.16) 

This remaining "sideband" term can be used to couple the atomic transition to the vibra- 
tional motion through a judicious choice of detuning, A = ±i^p, for some normal mode 



frequency Up. The first red sideband transition is obtained by setting A = z/^; that is, the 
laser is detuned one unit of vibrational energy helow (to the red of) the atomic transition. 
In this case, the resonant terms are 



■"> 



tT("i)/j.\ red Sideband, mode p ,„ Um , (m) , + ("^)^ / n a ^\ 

^i ^'^1 7. ; ' hVL^ri^j-^i^apaX + a],al ') . (3.17) 

(^ - "p) Up' 

This is of Jaynes-Cummings form |8l |20] , corresponding to excitation of the atomic transition 
with energy gap huA = f^h + ^t'p upon absorption of one laser photon at energy /icul, along 
with absorption of one vibrational phonon at energy hvp. Detuning the laser above (to the 
blue of) the atomic transition by one unit of vibrational energy (A = —Vp, for some mode p) 
generates the blue sideband transition, which has resonant terms 



hip) 

rr(m)/j.\ blue sideband, mode p ,„ Om , \ (m) , (m)x /o i o\ 

^i ^^) 7. ^ ' ^^o-n-^^WaX + Opai ') . (3.18) 



This corresponding to atomic excitation of the transition at energy ^a = ^l — ^^p upon 
absorption of one laser photon at energy fuvi,, along with emission of one vibrational phonon 
at energy hvp. 

The first red-sideband transition uses the ion itself as its own detector of motional quanta 
and this configuration comprises the first half of our two-stage detection model. As long as 
the coupling strength VLqT] is weak enough, no more than a single \g) -^ \e) transition will 
be excited for a given ion in the time during which the laser coupling is active. Although 
the approximations above simplify the interaction Hamiltonians and provide a means to 
understand the physical mechanisms at work, for finite detection times the full form of the 



interaction may be needed. As such, we will retain the "sideband" term of Eq. (3.15) as the 



interaction Hamiltonian, which has the form of a De Witt monopole coupling [21j, in order 



to allow for greater applicability of our results. In the examples given in Section |VI[ we will 
apply the long detection-time approximation at the end of the calculations. 

C. Excitation probabilities and correlation functions 

In the interaction picture, the full interaction Hamiltonian is 

Hj{t) := hQovY!<Pm{t)(^^'"\t) , (3.19) 

m 

where the prime on the sum indicates that only the ions being addressed by interaction 



lasers are included, and the unitless displacement operator (pmit) is defined in Eq. (3.7) 
The time evolution operator may be written formally as 



[/(T):=r|exp -^ /" dtHi{t) \ , (3.20) 



where T is the Dyson time-ordering symbol. The time-ordered exponential is defined by the 
time-ordering of its Taylor-series expansion, the first few terms of which are written here: 



U{T) = l + ^f dtHrit) + ^^j f dtij dhT{Hi{ti)Hi{t2)} + --- . (3.21) 



We will use these three terms in what follows. 

As discussed in the introduction, we wish to calculate correlation functions for simultane- 
ous detected excitations. Detecting the excitation of a given ion m corresponds to measuring 
the projector 

'Pm := \e)^Je\ , (3.22) 

where \e)^ is the excited electronic state of the mth ion. Vm acts trivially on all other 
electronic states and on all vibrational modes. In the experiment, the electronic state of 
the ion is determined using the technique of fluorescence on a cycling transition fS]. The 
excited state |e)^ is caused to make a dipole allowed transition to another auxiliary level 
which then decays back to the state |e)^ through spontaneous emission. If this transition 
is saturated, a very large fluorescent photon flux is easily detected. The measurement very 
nearly approaches a projection measurement of the operator Vm with an efficiency greater 
than 99%. 

These projectors commute for all m, so we can represent joint detection of the excitation 
of multiple ions as 

' mi---mM ' ' "11 ' ' ' ' niM ' yo.Zoj 

We will be interested here in simultaneous detection on at most two ions, 

'Pn.n = \e)^Je\ ® \eUe\ (3.24) 

(assuming m y^ n, since Vmm is just Vm), and initial states (in the interaction picture, at 
time t = 0) of the form 

Po = P®\9){9r , (3.25) 

where p is the vibrational state, and all of the ions are in the ground electronic state. The ions 
are assumed to be ordered in a linear array. Thus the subscripts mn are an implicit spatial 
index. This is indicated schematically in Figure [TJ The probability for measuring ion m in 
an excited electronic state after the interaction Hamiltonian is applied from time t = to 
time t = T is given by 

Pm := {U\T)VmU{T)) , (3.26) 

where the expectation value is taken with respect to po- Similarly, the probability that both 
ions m and n are in the excited electronic state after the evolution is 

Pmn := {U\T)VmnU{T)) . (3.27) 

Since the only outcomes of an actual measurement is a stochastic binary string indicating 
which ions lit up and which did not (ei, . . . , ctv), these probabilities correspond directly to 
classical expectation values over the entries Cm in this string: 

^ [CjT^J J^m ; 

Of course, other correlation functions (for three or more ions) may be defined analogously. 



^^ 
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FIG. 1: An illustration of the linear ion array with N ions. The index m runs from 1 to A^. This 
may be converted to a position label Xm as indicated, which labels the equilibrium position of each 
ion. The quantum degree of freedom, q^ describes small displacements from equilibrium. In (a) 
weak lasers couple the internal electronic state of the ion to the displacement from equilibrium 
of that ion. In (b) strong readout lasers drive fluorescence conditional on the electronic state of 
the ion. If an a given ion is in the excited state, it fluoresces (giving the result 1) or it does not 
fluoresce (giving the result 0). 

IV. EXCITATION PROBABILITY CALCULATIONS FOR GENERAL STATES 

We proceed to calculate Pm and Pmn to second order in the Dyson series expansion, 



Eq. (3.21). The coupling lasers are assumed to be weak enough to justify this perturbative 
treatment. The zeroth-order approximation Pm vanishes since the initial state has no ions 
excited. The first-order term is 



p(i) 



1 



T 



dh / dt2 {Hi{ti)VmHi{t2))tom 



^0 

T 



(Qovf dh dt2jg\a^:-\h)\e)^Je\a^:-\h)\g)^{Mti)Mt2)), (4.1) 



resulting in 



P« = (fior/)'/ dh dt2e-'^^'^-''\4>Uti)4>m{t2)) 



(4.2) 



where expectation values are now over the vibrational modes only unless otherwise indicated. 



,(2) 



(m) (m) I 



The second-order term P^' = because J[g\Hj{ti)Hi{t2)\e)^ ~ m(5'l'^± "^i f)^ = — i-^-, 
no two applications of the electronic raising/lowering operators can lower the excited state 
to the ground state. 



We're now prepared to tackle the spatial correlation function (3.27). Once again, the 
zeroth-order term Pmn = because no ions are electronically excited to begin with. Also, 
considering that Vmn projects onto the excited electronic states of two distinct ions while 
each Hi{t) can only raise one ion at a time, ^(^'1 ®n(fi'l-^^('^i)k)m® k)n = 0; making Pml = 0, 
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as well. Therefore, we must go to second order: 

cT pT rT pT 



P!nl = ^,f dhj dtj dtj dU{T{Hj{h)Hjih)}V^a{Hj{ts)Hj{t,)}),,,^, 

X T{H'r\h)H'f\u) + Hf\t,)H^r\u)}U.i . (4.3) 



4/i^ ./o 



The antitime-ordering symbol T acts like T but instead orders the terms from earliest to 
latest. Trading integration variables (ti -v^ t2 and/or t^ <^ t4) and minding the (anti)time- 
ordering, we can simplify this to 

Pl^l = ^ ^^rf^t(T{/jf )(t,)/ff )(tl)}P^nT{/fr (t3)iff)(t4)})total • (4.4) 

The terms from the Hamiltonians that survive the projection and contraction with the 
electronic ground state are 

f {i/f )(t2)i^f ^(ti)} =^ {hQoV?T{Mt2)Mti)}e~''''^e-'^'^a'r^a^^^ (4.5) 
and T{i/f )(t3)//f ^(t4)} =^ (;ifior/)'T{0„(t3)0„(t4)}e+^^*^e+^^*vf Vf . (4.6) 
Plugging these in gives 



Pj^l = (Qov)' [ dHe-'^^'^+'--'^~'^\f{<PUt2)Uti)}T{Mt3)MU)}) ■ (4.7) 



Several ways exist for expressing the four-point correlation function in Eq. (4.7) in a con- 
venient form by use of Wick's theorem. These are included in the Appendix. One form is 
particularly useful, though, and that it applies when the state has a Wigner function that is 
a Gaussian with zero mean (i.e., a "Gaussian state"). In this case, we may write (repeated 



from Eq. (A. 18)) 



(T{0^(t2)0n(tl)}T{0„(t3)0„(t4)}) '""""""") (0m(^2)0m(i3)) (0n(il)0n(i4)) 

state 

+ {(pm{h)MU)){Mtl)4>m{h)) + {r{(t)m{t2)Mtl)}){T{(PUt3)MU)}) ■ (4.8) 



Comparing this with Eq. (|4.2), we see that the first term will always give simply PmPn, and 



the second will always give a term similar to this, but with a different geometric factor. Thus, 



once we have Pm, we need only ever explicitly calculate the last two terms from Eq. (4.8) to 

§6t i-mn- 

A. Long-time interaction 

Before moving on, let's have a look at what affect a long-time interaction would have 
on the evolution. In this case, we can employ the rotating wave approximation and use 



11 



Eq. (3.17) as the interaction Haniiltonian. In this case, Eq. (3.20) can be evaluated to 

UiT) 



RWA 



(A = i/p) 



exp 



zT{nov)J2 



'h 



(p) 



(m) , t (»^)^ 



(4.9) 



There is a balance to be maintained, here, though. On the one hand, the detection time T 
needs to be long enough so that the rotating wave approximation is valid. This is the 
requirement that T ^ i>~^ (and we assume that A ~ i/). On the other hand, it needs to be 
short enough so that we can use perturbation theory, which requires T <C {Qori)~^. When 
both of these conditions are satisfied simultaneously, that is. 



i^"' < T < (noriy' 



we can expand the exponential in Eq. (|4.9|), as in Eq. (|3.21|). This gives 

rp2 



p(i) 



RWA 



{HiV^Hj) 



total 



Jp)2 



(4.10) 



(4.11) 



which agrees with what is obtained from selecting only the resonant terms directly from 



Eq. (4.2). A similar result holds for the second-order function: 

rp4 



p(2) 



RWA 



(A = !/p) Ah^ 



:^/^™n^/>totai = T\^o^y 



(fe^^fei^^)^ 



f^p 



lp0'p(^p^ 



aLdLdndp) 



(4.12) 



If we were to allow for different detunings on each ion, A^ 
could be generalized: 



'^p; ^n 



Upi, then this result 



p(2) 



RWA 



(A„ 



,.p, An = iyp>) Ah* 



> — (H]V„,nHf) 



total 



T^Qov)' 



(Mdp')\2 
yUm Un ) 

^fipflpi 



t t 
Ojpdpidp'Ctp 



(4.13) 



This interaction would require that a second laser beam be detuned from the first, an 
experimental complication we won't consider further. 

For typical values of the parameters (z/ ~ 1 MHz, Qq ~ 100 kHz, and u ~ 0.01), Condi- 
tion (4.10) requires that 10~^ s <^ T <^ 10^'^ s, a rather narrow band in which to operate 



well within both regimes. Given this limitation, in the next section, when we calculate Pm 
and Pmn for any Gaussian state as a function of its covariance matrix, we will retain the 
perturbative condition T <C {Qori)~^, but we will not use the rotating wave approximation, 
in order to allow for a more general class of interactions, including ones for which T < i/~^. 
We will also assume that both interaction lasers have the same detuning A for experimental 
simplicity. 



V. EVALUATION FOR GAUSSIAN STATES 

A. Two-point functions 

As discussed in the Appendix, when the vibrational state is a zero-mean Gaussian, all 
measured correlation functions can be evaluated from the two-point functions {(f)m{t)4'n{t')) 
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and (T{0m(t)0„(t')}). Let's define the two-time-dependent matrices 

T(t,t'):=(0(t)0(t')^> (5.1) 

and f{t,t') := {T{(t){t)(l){t'f}) (5.2) 

(5.3) 

for which T„„(t,t') = {4>mit)4>nit')) and T„,„(t,t') = (T{0m(t)0„(t')}). It's easy to see from 
the entries that 

^<-')^ {;!:;::!• ;;::::; 

since (f)m{t) is Hermitian. Thus, we really only need to worry about T(t, t') for the moment. 



Using Eqs. (3.7) and (3.2), we can transform to the normal-mode two point function instead: 



T(t, t') = B^A-i/^ ( [a(t) + a(t)] [a(t') + a{t')f) A^^/^B , (5.5) 



where a(t) := (ai(t), . . . ,ai^{t))'^ is a column vector of interaction-picture lowering opera- 
tors for the normal modes, a(t) := {a[{t), . . . ,a'j^{t))'^ is the equivalent column vector of 



interaction-picture raising operators, and B and A are defined through Eq. (3.2). We can 
write the time-dependence of these vectors in a very compact form: 

a(t) = E(t) o a , (5.6) 

where 

E{t):={e-"''\...,e-"'^Y (5.7) 

is a column vector of time-dependent coefficients, and the symbol o represents the Hadamard 
product (element- wise multiplication). Similarly, 

a(t) = E(-t) o a . (5.8) 

Any Gaussian state with zero mean is uniquely defined by its covariance matrix. While 
many varieties of covariance matrix can be defined [19j, an obvious choice here would be 
to use the two matrices (aa-^) and ^aa-^). We can get the other combinations by noting 
that (aa ) = (aa^) , and (aa ) = (aa^^) + 1. Thus, we can go one step further and define 
a matrix of coefficients 

E(t,t'):=E(t)E(t')^ (5.9) 

such that Ers{t,t') = ^-'^i^rt+vst )_ Using this shorthand, we can isolate the time dependence 



from the expectation value in Eq. (5.5) and write the result in terms of the initial covariance 
matrix: 



K(t, t') := ^ [a(t) + a(t)] [a(t') + a(t') 

= (aa^) o E(t,t') + (aa^) o E(t, -t') + (aa'^> o E(-t,t') + (aa^) o E(-t, -t') 
= 2 Re [(aa^) o E(t, t')] + [(aa^) + l] o E(t, -t') + (aa^) o E(-t, t') . (5.10) 
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Considering Eq. (5.4), we're going to need the complex conjugate of Eq. (5.10). Element 



by-element evaluation will reveal the following equivalences: 

[(aa^> o E(t, -t')]* = Kaa^> + 1] ° E(-t,t') , 
[(aa^>oE(-t,t')]*=Kaa^>-l]oE(t,-t'), 

leading to 

[(aa^) o E(t, -t') + (aa'^) o E{-t,t')]* 

= (aa^) o E(t, -t') + (aa^) o E(-t, t') + 1 o [E(-t, t') - E(t, -t')] 
= [(aa^) + 1] o E(-t,t') + (aa^) o E(t, -t') . 



(5.11) 
(5.12) 



(5.13) 



This is the same as the last two terms in Eq. (5.10), except that the E- matrices have been 
exchanged. Therefore, we can define 



E{t,t'):- 



E(t, -f) if t > t' , 
E(-Lt') ift<t'. 



(5.14) 



Notice the subtle difference between T(t, t') and E(t, t'). The ring-notation is consistent with 
its purpose — to describe a function associated with time-ordering — but not necessarily with 
the details of the definition. Recalling Eq. ( |5.4 ), the time-ordered expectation value can 
now be written succinctly: 



K(t, t') := (t { [a(t) + a(t)] [a(t') + mf}] 

= 2Re [(aa^) o E{t,t')] + [(aa^) + l] o E(t,t') - 

Pluggin into Eqs. (5.5) and ( |5.4 ) gives 

T(t,t') = B^A-^/^K(t,t')A-'/'B 
t(t,t') = B^A~^/^K(t,t')A-'/'B 



(aa^)oE(t,t') 



(5.15) 

(5.16) 
(5.17) 



with K(t,t') defined in Eq. (5.10) and K(t,t') in Eq. (5.15). 



B. Probabilities in terms of the covariance matrix 

We have successfully consolidated all of the time-dependence into E(t, t') and E(t, t'). This 
makes evaluation of Eqs. (4.2) and (4.7) dependent only on integrals involving elements of 



these matrices. The following integral will be useful: 



1 
T 



T/2 



dte 



iujt 



smc 



T/2 



ujT 



(5.18) 



where 



smc a; := 



X ^ sin X if x 7^ 
1 if a; = 



(5.19) 
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and we have symmetrized the hmits of integration by setting the laboratory clock 
appropriately — a passive operation that does not affect the physics of the experiment. 



Eq. (5.18) is a bandwidth-limited Fourier transform. In anticipation of future calculations, 



let's define 



f.T/2 i'T/2 



1 /•-'/^ /-J/^ 

Srs{uj, J) :=— dt dt' e'^^'^'^'^'^Ersit, t') 

1 J-T/2 J-T/2 



smc 



-T/2 J-T/2 



Sine 



(5.20) 



We can now collect the Sr^ioj.oj') elements into a matrix Siiuj.uj'). We can also define 



SrAuJ.uj') :-- 



1 j-T/2 I'T/2 

2^2 



/ dt dt'e'^^'+^'^Ers{tX) 

J-T/2 J-T/2 



(5.21) 



and a corresponding matrix S{uj,uj'), although we will leave it unevaluated for now. 

m ■ 



Using these tools, let's calculate Pn 

rT/2 
/ rft2e-^^(*^-*^)T„„(ti,t2 

J-T/2 



r.T/2 pT/2 

Pi^) = (novf I dh 

-T/2 J-T/2 



(fior/)^e^B^A-V4 



rT/2 nT/2 

/ dti / cit2e-*^(*^"*^)K(ti,t2 

J-T/2 J-T/2 



A-i/^Be. 



(aa'^) o S(-A, A) + (aa^) o S(A, -A) 



[(aa^) + 1] o S(-A, -A) + (aa^) o S(A, A) 



A-i/^Be. 



(5.22) 



where e^ is a unit column vector used (twice) to pick out the correct element from the 



matrix. We should also point out that Eq. (5.22) also holds for non-Gaussian states, for 



which a covariance matrix can also be defined even though it does not specify the state 
completely. 

The correlation function P„m can be calculated similarly. In this case, it is important 
that the state be (zero-mean) Gaussian because we will use the simplification provided by 



Eq. (4.8). As already discussed. 



(Term 1) = P^nPn 



(5.23) 



The second is almost the same, except it involves a different element of T(t, t'). An analogous 
calculation to the one above shows that 



(Term 2) = T^(fio^)' 



£b^a-i/^ 



(aa^) o S(-A, A) + (aa^) o S(A, -A) 



[(aa^) + 1] o S(-A, -A) + (aa^) o S(A, A) 



A-i/^Be. 



, (5.24) 
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the only difference (besides tlie squaring) being tfie presence of e„ at tlie end, instead of e^. 
The third term is more complicated, due to the time ordering. Nevertheless, it can be 
written 



(Term 3) = (fio^/)^ 



T/2 /•T/2 

dti I 

T/2 J-T/2 



/.T/2 /.T/2 

/ dh / rft2e*^(*i+*^)T„„(ti,t2 

J-T/2 J-T/2 

/.T/2 /.T/2 

/ dti / rft2e*^(*i+*^)K(ti,t2 

J-T/2 J-T/2 



A-^/^Be„ 



elB^A-i/^ 



(aa^) o S(A, A) + (aa^) o S(-A, -A) 



+ [(aa^) + 1] o S(A, A) + (aa^) o S(- A, -A)* 



A"i/%e. 



(5.25) 



Evaluating this term boils down to evaluating Eq. (5.21). The sum of Terms 1, 2, and 3 

(2) 

gives Pmn for a zero-mean Gaussian state. 



VI. EXAMPLES 

In order to demonstrate the usefulness of the correlation function, we'll compare a thermal 
state with a given average phonon number to a corresponding uniformly squeezed state with 
the same average phonon number. We will also assume that the interaction is active on a 
timescale long compared to the period of vibrations of the ions — that is, T S> z/~^. In both 
cases, the probability of excitation for any given ion is the same, but, as we shall see, the 
correlations are stronger in the squeezed state versus the corresponding thermal state. Our 
measure of correlations will be 



fn 






(6.1) 



which is a normalized correlation function akin to those defined in Eqs. (2.5) but is based 



on detection probabilities, rather than underlying quantum correlations. The connection to 
the quantum state, of course, was made in the preceding sections. 



A. Thermal state 



A thermal state at temperature r is a zero-mean Gaussian state in which 



{aras) = , 

(alas) = UrSrs , 



where 



Ur '■= 



1 

pPhVr ][ 



(6.2) 
(6.3) 



(6.4) 
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is the average number of phonons in mode r, and /3 = (/cbt) ^- Plugging into Eq. (5.22) 
gives 

Pii) = T2(fior/)'5^^{(np + l)sinc2[(A + z/,)f]+n,sinc2[(A-z/,)f]|. (6.5) 



The sinc^-function is sharply peaked at A = ±i/p, and in the long-time limit {T ~^ v ^), we 
have 



-1 Ti'-P-'^ 



^(rt2 

where the Kronecker-5 symbol satisfies 6q = 1 and 6^ = for x 7^ to within a band- 
width of approximately T~^. This is gives the same results as would be obtained after using 



the rotating wave approximation, as described in Section IV A Choosing a particular nor- 
mal mode frequency as the detuning (A = Up) accords with the result calculated directly 
from Eq. (4.11). In this case, the detection probability has a geometric factor yUp ' bm\ 



corresponding to the position of the ion within the normal mode p being addressed and is 
proportional to the average number of phonons Up in the mode for red sideband detuning 
(and np + 1 for the blue sideband), as expected. 

Moving on to the correlation probability Pmn, we can evaluate Eq. (5.24) to 






(Term 2) = 
In the long interaction time-limit, this term behaves similarly to P, 



-\ 2 



^^!^{(np + l)sinc2[(A + z.p)f]+npSinc2[(A-^p)f]| 
j^p y ) 



■ (6.7) 



(1). 

m ■ 



(Term 2) il^^ T\^Q^f 



E 



l(p)l(p) 



{up + l)(27r)^(5(A+j,p) ■\- np{2ny6(^A~up) 



(6.8) 



In order to evaluate Eq. (5.25), we need to evaluate S(±A, ±A) from Eq. (5.21). Using 



Eq. (6.3), E(ti,t2) is diagonal, meaning we only need to calculate 

f.T/2 

/ dt. I 

J^2 



1 rT/2 pi/z 

Spp{±A, ±A) = — / dh / dt2 e±*^(*^+*^)Ppp(ti, t 

J-T/2 



-T/2 

1 rT/2 

2^2 



T/2 

rT/2 

dtl / (lf^Q±iMtl+t2)^-iUp\tl-t2\ 

-T/2 J-T/2 



O rT/2 ft 

T J-T/2 J-T/2 

Since T ^ z/~^ ~ A, we can change the integration limits with T -^ 00: 



^pp(±A,±A) 



2^2 



dti / dt2e^^^(*^+*^)e-*''^(*i-*^) 



(6.9) 



(6.10) 
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With the integration now over the entire half-plane defined by t2 < ti, we can rotate our 
integration axes using 

ti = U + V U = -{ti + t2) 

t2 = U-V V = -{ti - t2) 



(6.11) 



to obtain 



/■oo noo 

Spp{±A, ±A) = 2 / dv due^'^^^'e-'^'"'''' ~ 5a5^, ^ . (6.12) 

Jo J~oo 



Thus, in the thermal case, 



(Term 3) il^iL^l, q . (6.13) 



Consolidating these results, we have 



Pr 



(2) 



•^'^"- p(l)p(l) ~^' ^^"^^^ 

for any detuning A = ^Up, any nonzero temperature, and any choice of ions {m, n) to 
measure. This is consistent with the expected results for a second-order correlation function 
for a thermal state 1191 . 



B. Uniformly squeezed normal modes 

The state to be considered in this section is uniformly squeezed in all normal modes. 
Such a state has 

{a^as) = -^/n{n + 1) 6rs , (6.15) 

(alas) =n6rs , (6.16) 

where n = sinh r is the mean phonon number of each mode (assumed the same for each 
mode) as a function of the squeezing parameter r. 

Proceeding as for the thermal state, we have nearly the same expression for Pm ■ 

N ,(P)2 

P}n^ = T\noVr E ^ I (^ + 1) ^i^^' [(^ + ^r>) ^]+n sinc^ [(A - z.,) f ] 

p=i V^^p 

- 2^n{n + 1) sine [(A - z/p)f ] sine [(A + Vp)\\ \ . (6.17) 
In the long detection time-limit, however, the last term makes no difference, and 

(1) (^»^~^), p(l) ,. .ON 

m, squeezed jn, thermal ' V / 



>(1) 



where -Pi thermal is -^1- (6.5) for a thermal state with temperature chosen to make n = np in 



Eq. (6.4) for a chosen detuning of A = //„ 
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(2) 

An analogous calculation for the second term of Pmn shows that 



(Term 2) = T^VLqT])' 



hip) hip) 



Y^ ^!!i^ { (n + 1) sinc^ [(A + z/^) f ] + n sinc^ [(A 



^p)f] 



-\ 2 



2v/n(n + l)sinc[(A - i/p)|] sine [(A + i/p)|] 



. (6.19) 



In the limit of T ^ z/ , this gives no change from the same term in the equivalent thermal 
case: 

(Term 2),„^, ^^^^^ (Term 2\^^^^^, , (6.20) 

In the thermal case, the third term vanished in the long detection time-limit. In the squeezed 
case, it does not, and this generates the difference in the correlation functions: 



{r»i.-i) ^4 4 



(Term 3) ^ "^ ^ T^inovY 



u('')h(^) 






(6.21) 



This shows that even though local measurements have the same excitation statistics as the 
equivalent thermal state, simultaneous measurements of two ions do not — the correlations 
are stronger in the squeezed case, as should be expected. For the uniformly squeezed state, 
we have 

p(2) . 

/.„-^;(^ = 3 + -, (6.22) 

for any detuning A = it'p and any choice of ions (m, n) to measure. The value depends on 
the squeezing parameter through n. The reason the two do not agree in the limit n — > 
is that in that case Pm = Pn = 0, so fmn is undefined. The behavior of this function, too, 
is consistent with the expected results for a second-order correlation function for such a 
squeezed state [19j. 



VII. DISCUSSION AND CONCLUSION 

There are really two pieces of information that can be gleaned from the correlation func- 
tion Pmn- The first is the structure of the normal modes. This structure can be probed by 
detuning the interaction laser to the desired mode's resonant frequency i>p and reading out 
ion m and n. Doing this for either the thermal state or the uniformly squeezed state from 
the previous section give correlations as in Figure |2} 

The other piece of information obtainable from Pmn is its relation to the single-ion detec- 
tion probabilities Pm and Pn- This information is embodied in the normalized function fmn, 
and can be used to probe more general characteristics about the state. For instance, in 
the comparison of squeezed and equivalent thermal states, with equivalence defined as an 
equal average number of phonons in the mode being detected (i.e., n = n^ for A = Up), 
many of the parameters cease to be important (detection time, coupling constant, geometric 
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FIG. 2: Correlation functions Pmn for the normal modes of 10 ions in the case of the examples 
considered in Section VI The top row represents Pmn for a detuning corresponding to normal 
modes 1 through 5; the bottom row, modes 6 through 10. Since Pmn oc PmPn for both types 
of states, except in the case of the center-of-mass mode (mode 1), which is uniform anyway, the 
colors are normalized such that white represents the maximum value, and black represents the 
minimum in each case. Thermal and uniformly squeezed states both give the same results under 
this normalization. The diagrams show two things: (1) that the structure of the modes revealed by 
the correlations, and (2) which pairs of ions (white boxes) are most useful for measurement when 
detuning to a given mode frequency. 

terms, etc.), and information about the basic nature of the state — in this case, thermal or 
squeezed — is revealed. 

There is still much work to be done in exploring properties of correlations in a string 
of trapped ions. The work presented here is designed to be a significant start in that 
direction. We began with a definition of a measurable correlation function and single-ion 
detection probabilities and connected these to two- and four-point functions of a general 
quantum state, all in the perturbative regime. Next, we specialized to Gaussian states and 
provided explicit formulas for these probabilities in terms of the state's covariance matrix 
elements. The examples of thermal and squeezed states were compared, and a normalized 
correlation signature contrasted in each case, which shows agreement with standard results 
for correlation functions used in quantum optics. Open problems include vast opportunities 
to generalize these results to other interesting states, such number states, coherent states, 
superpositions versus mixed states, and many others. The structure of a string of trapped 
ions approximates a scalar field, but the deviations from this approximation could be a 
source of interest, as well. Finally, a completely unexplored avenue would be to look at the 
behavior of these correlation functions as the trap frequency is modulated, to see whether 
such a modulation has a detectable correlation signature. 
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APPENDIX: WICK'S THEOREM 

Wick's theorem theorem is often stated as a result for time-ordered expectation values of 
the ground state of a multimode bosonic or fermionic system, as in quantum field theory |22] , 
but it can actually be used with any state and any prescribed ordering of the operators. 
Wick's theorem is commonly stated as [22]: 

T{(f)mAti) ■ ■ ■ 0m^^(iAf)} = :0mi(ii) " " " 0mjv(^jv): + :(all possiblc contractions): , (A.l) 

where the colons in : (operators): place the contained operators in normal order (i.e., all 
raising operators are to the left of all lowering operators), a "contraction" of two operators 
is written as (and defined by) 



[ijrn(tj),ijl(tk)] if tj > tk 






where ipmit) and iplnit) are defined in Eq. (3.9), and "all possible contractions" means every 



unique way of contracting pairs of operators together in this fashion. 

We can make two generalizations: to anti-time-ordered products and to products with 
no time ordering. Since normal ordering already prescribes and order for the uncontracted 
operators, the only place time ordering shows up at all is in the contracted terms. As such, 
we define a second type of contraction for this purpose: 

, (. ^ , (. ^ ._ j[Mtj),i^i{tk)] iitj <tk, , . 

I 1 ( ['4^n{tk),Vm{tj)] if tj > tk ■ 

Finally, if the operator order is not prescribed by either type of time ordering (and is to be 
taken as given), then we shall define 



^m{tj)Mtk) ■■= li^mitj), ijiitk)] . (A.4) 



Generalizing the usual inductive method of Ref. [22], it's straighforward to generalize Wick's 
theorem to the following: 

^{0mi(i^i) ■ ■ ■0mAr(i7v)} = ■(pmAh) " ' ■0mjv(%): + :(all posslblc coutractious) : , (A.5) 

where V{- ■ ■ } prescribes a time ordering for (some of) the operators within it, and the 
contraction terms now must respect this ordering, using one of the above definitions in each 
case according to the prescribed ordering of the two operators involved. 



For our purposes, we are interested in the four-point function in Eq. (4.7). This also 



serves as a particularly good but simple example of how to use Eq. (A.5). We'll use the 
following abbreviations to save space: 

Mtl) -^ 0' , (Pmih) -^ 0' , 0m(t3) ^ 0' , MU) -^ 0' ■ (A.6) 
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We then can use Eq. (A. 5) to make the following expansion: 



+ :^ 



I I 



: + :02cF0^0^: + :<P^0^: 
I 1 ^ ^ I 



r 



n, 



+ :<p0Vj.^: + :02<FJ)30^: + :<^2<^i<^3<^^: 
I I I I I I I I 



(A.7) 



Notice how the contraction type is determined by whether the terms are prescribed to be in 
time order {(fr" and 0^), anti-time order (0^ and 0^), or as written (all other combinations). 
We need the expectation value of each of these terms. Since contractions of position operators 
are c-numbers, they can be taken out of the normal ordering and out of all expectation values, 
giving 

(T{020i}T{030n) = (:</>'0V'0':) 

+ (Fi''(:0V':) + {■■<f<f:)W + WW 
I I I II ^ I I 

+ W{-4'<P':) + {■4^<P':)W + WW 
I I I II ^ I I 

+ <P^\-4^':) + {■4'<P'-)W + 0V'<F0' • (A.8) 



The grouping of the terms in Eq. (A.8) suggests an alternate way of writing the expectation 



value of this expression. Using Eq. (A. 5), we have 

(0V)(0V') = (:0V:)(:0V':) + W{-4'<P'--) + {^'^t^'-^W + WW , (A.9) 



a2 j4\ / j1 j3\ 



/,2'i4 lili/S 



0^)(0V^) = (:0V:)(:0V^:) + <^T(:0V^:) + (:0>^:)<^V^ + r<^^<^^<^' 



a3^4. 



2^1/.^3 i,4.\ 



/,2 J,l.\ ^3 114 



l2i1}i3[i4 



(r{0^0^})(T{0^0n) = {■■<i>V:){-4'r--) + rri-.rr--) + {--vvwr + rvr^ 



Using this, we can write 



(A.IO) 



(A.ll) 



(T{020i}T{030n) = (0'0')(0V') + (0'0')(0V') + (T{0Vi})(T{030n) 

+ (:0VW:) - (:0V:)(:0V':) - (:0V:)(:0V':) - (:0V:)(:0V:) • (A.12) 



Recalling the abbreviations (A. 6), either Eq. (A.8) or Eq. (A.12) may be used to calculate 
Pm-n to lowest nontrivial order. 



Gaussian States 

If the state has a Wigner function that is a Gaussian with zero mean, then all of the 



normal-ordered terms in Eq. (A.12) must cancel out in order to agree with the result known 



as the "generalized Wick's theorem" from Ref. [23]. As it is stated in Ref. [23j, the theorem 
applies only to ground-state expectation values of raising and lowering operators, but it can 
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be generalized to a much larger class of expectation values. Let Kj G {op, aj, \ p = 1, . . . , N} 
be any raising or lowering operator. The generalized Wick's theorem states that any ground- 
state expectation value of an even number of such operators can be written in the following 
form: 

(0| fi:ifi:2fi:3 ■ ■ ■ fta™ |0) = ^ (0| K1K2 |0) (0| /tgM |0) ■ ■ ■ (0| fi:2n-iK2n |0) , (A.13) 

Pd 

where the sum is over all distinct permutations Pd of the 2n indices that preserve the operator 
ordering within each pairing on the right — i.e., all permutations which give a distinct product 
of expect at ion- value pairs {niUm) and for which / < m in each of these pairs. Also, any linear 
functions Fi = ^^ ^iji^j of the raising and lowering operators will separate in a similar 
fashion: 

(0| F1F2F3 ■ ■ ■ F2„ |0) = 5^ (0| F1F2 |0) (0| F3F4 |0) • ■ ■ (0| F2n-lF2n |0) , (A.14) 

Pd 

as can be verified by direct calculation. Furthermore, it turns out that this property holds 
for all zero-mean Gaussian states, both pure and mixed. To show this, one notes that any 
such Gaussian state can be written as the partial trace of a pure Gaussian state of twice 
as many modes [2^: p = tr^ |x)(xl- Such a Gaussian pure state |x) is unitarily related 
to the vacuum state on the doubled mode set — i.e., \x) = U^ |0,0). However, this unitary 
transformation can be viewed in the Heisenberg picture as a symplectic linear transformation 
on the (extended) vector of canonical operators [25], viz. U^HjU^ = Ylk-^^Vk ='■ Gj, where 
Kj is any raising or lowering operator in the original set, while ifk is any such operator in 
the extended set. From this, we can write expectation values as 

{ki--- Kn) = tr(pKi ■■■/€„) 

= tlA[trB{\x){x\)filf^2fi3 ■ ■ ■ /«2n] 

= (0,0|f/t«:i[/^...f/t«:„?7jO,0) 

= (0,0|Gi---G„|0,0) . (A.15) 

Thus, we have 

{KIK2K3 ■ ■ ■ K2n) = (0, 0| G1G2G3 ■ ■ ■ G2n |0, 0) 

= (0, 0| G1G2 10, 0) (0, 0| G3G4 |0, 0) ■ • ■ (0, 0| G2„-iG2„ |0, 0) 

= {k,iK2) {K3K4) ■ ■ ■ {K2n-lK.2n) , (A. 16) 



where we applied Eq. (A.14) in the middle step. Thus, the generalized Wick's theorem holds 
for all Gaussian states with zero mean. Furthermore, we may take linear combinations of 
Kj operators, and the theorem still holds: 

(F1F2F3 ■ ■ ■ F2n) = 5^(FiF2)(F3F4) ■ ■ • (F2„_iF2„) . (A.17) 



If the state in Eq. (A. 12) is a zero-mean Gaussian, then Eq. (A.17) applies, resulting in only 



the first three terms on the right being kept: 

(f{0„(t2)0n(tl)}T{0„(t3)0„(t4)}) '"'""''"'> {(pm{t2)(pm{t3)){<Pn{tl)(pn{h)) 

state 

+ (0m(t2)0n(t4))(0n(tl)0™(t3)) + (T{0„(t2)0n(tl)}) (T{0^(t3)0n(t4)}) • (A.l^ 
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Comparing this with Eq. (4.2), we see that the first term will always give simply PmPn, 



and the second will always give a term similar to this, but with a different geometric factor. 



Thus, we need only ever explicitly calculate the last two terms from Eq. (A. 18). 

We also note that if we restrict to Gaussian vibrational states, the generalized Wick's 
theorem says that all higher-order joint probabilities (those with more excited ions, like Pimn, 
and/or higher-order calculations) will all decompose into integrals and sums of expectation 
values of the form (0m(^i)0n(^2)) and (T{0m,(^i)0n(^2)}) for arbitrary m and n. (Note that 
the antitime-ordered case is just the complex conjugate of the time-ordered case). When 
restricting to Gaussian states, then, only these two types of expectation values are required. 
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